GPU-accelerated locally injective shape deformation

ABSTRACT

A computer-implemented method for shape deformation, comprising: obtaining an image; obtaining one or more positional constraints; determining a shape deformation/image map, based on the positional constraints; applying the shape deformation to the image; and outputting the image.

The present invention relates to a method and a device for deformingshapes in images.

INTRODUCTION

Shape deformation is a classical problem in computer graphics andanimation that comes in a variety of forms and settings. Essential toall variants is the requirement to compute a map between differentspaces with some desirable properties. Common choices include mapsbetween curved manifolds embedded in R³, volumetric maps between solidshapes in R³, and mixed dimensional maps where, for example, a curved2-manifold embedded in R³ is mapped to R² (a parameterization).

Depending on the precise user requirements, even a simplified settingusing planar maps, where both the source and target domains are subsetsof R², is challenging to deal with.

Some requirements such as smoothness are easily achievable while other,such as local injectivity or strict bound on the induced metricdistortion are hard, sometimes impossible to guarantee. The goal is tobe able to compute a map that adheres to these positional constraintsyet at the same time is visually plausible. Plausibility is usuallyquantified by smoothness and the amount of metric distortion, which canbe either minimized on the average, or bounded above by someuser-defined amount. Moreover, the map is expected to be locallyinjective, preventing degeneracies of the image and local overlaps whileallowing global ones. These requirements pose a great computationalchallenge since they usually lead to a nonconvex optimization problemwith many degrees of freedom for which finding an optimal solution (oreven just a feasible one) is often intractable.

Recent advances in the computation of locally injective and/or boundeddistortion maps allow artists to produce high quality results withgreater than ever success rates but the performance of current methodsis still significantly worse compared to simple linear methods (WeberOfir, Ben-Chen Mirela, and Gotsman Craig. 2009. Complex BarycentricCoordinates with Applications to Planar Shape Deformation. ComputerGraphics Forum 28, 2 (2009), 587-597).

It is therefore an object of the invention, to provide a computationallyefficient method and a device for shape deformation of an image thatproduces visually plausible results.

This object is achieved by a method and a computer program productaccording to the independent claims. Advantageous embodiments aredefined in the dependent claims.

According to the invention, the deformation problem is cast as anonlinear nonconvex unconstrained minimization problem that is performedin a reduced deformable subspace of harmonic maps. The user is allowedto pick among a variety of smooth energies that measure isometricdistortion.

The method according to the invention broadens the scope of previousresearch on harmonic maps from simply-connected to multiply-connecteddomains (disks with finite number of holes) such as those appearing inFIGS. 2, 4, 9. In particular, the method extends the harmonic subspaceof the Bounded Distortion Harmonic Mappings (BDHM) method ([Chen andWeber 2015] Chen Renjie and Weber Ofir. 2015. Bounded distortionharmonic mappings in the plane. ACM Transactions on Graphics 34, 4,Article 73 (2015)) to multiply-connected domains, as well as generalizesthe bounded distortion theorem of BDHM. Furthermore, a much tighteranalysis of the Lipschitz constants provided in BDHM is provided,leading to more accurate validation of the map and faster convergence.

A custom-made Newton solver is constructed with features efficientlyaddressing the deformation problem. The quality of the produced maps issuperior and comparable to the other methods that operate within theharmonic subspace. However, the inventive method is orders of magnitudefaster than existing techniques. The optimization is based on asecond-order Newton approach which ensures that the convergence rate ishigh. In general, a single Newton iteration can be slow due to variousreasons such as the need to approximate the Hessian numerically usingfinite differences (when analytical expressions are unavailable), aswell as the need to modify the Hessian such that it becomes positivedefinite, which is necessary to ensure convergence to a local minimum.

The inventive solver overcomes these issues due to several uniqueproperties, such as having closed-form expressions for the map, and itsdifferentials, and an analytic modification scheme for the Hessian toensure its positive definiteness. Combining these benefits leads to aNewton iteration that is considerably faster than a standard Newtoniteration. Together with the small number of iterations that aretypically required for convergence, the total running time of our solveris significantly shorter compared to state-of-the-art solvers. Moreover,every iteration of our solver is embarrassingly parallelized andperfectly suites implementation on massively parallel graphicsprocessing units due to the structure and regularity of the problem.This allows for unprecedented performance.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawing(s) will be provided by the Office upon request and paymentof the necessary fee.

These and other aspects and advantages of the present invention willbecome more apparent when considering the following detailed descriptionof an embodiment of the invention, in connection with the drawing, inwhich:

FIG. 1 shows a schematic flowchart of a method for image deformationaccording to an embodiment of the invention.

FIG. 2: shows a schematic drawing of a knot.

FIG. 3: shows notations.

FIG. 4: shows harmonic deformations of a triply-connected domain.

FIG. 5: is a table (Table 1) listing several possible choices for smoothisometric energies in different embodiments of the invention.

FIG. 6: shows a comparison of isometric energies E_(iso) (left) andE_(exp) (right) used according to different embodiments of theinvention.

FIG. 7: shows a comparison of Lipschitz constants according to anembodiment of the invention and according to prior art approaches.

FIG. 8: shows, in pseudo-code, an algorithm (Algorithm 1) for a Newtonsolver used to determine an image deformation map according to anembodiment of the invention.

FIG. 9: shows locally injective isometric deformation of Rex, adoubly-connected planar domain (left).

FIG. 10 shows a comparison with BDHM using the same harmonic subspace.

FIG. 11 shows the effect of changing the number of Hessian samples |H|on the algorithm according to an embodiment of the invention, whenapplied to the depicted deformation of the deer.

FIG. 12 shows how the runtime of the inventive GPU-based modified-Newtonsolver is spread across the 4 main stages that compose a singleiteration.

FIG. 13 (Table 2) shows a comparison of runtime (ms) of a singleiteration.

FIGS. 14 and 15 show a visual comparison of our algorithm againststate-of the-art methods

FIG. 16 shows additional results minimizing E_(iso) ^(ƒ) rendered athigh-resolution.

DETAILED DESCRIPTION

According to an embodiment of the invention, an interactive deformationmethod is described that is driven by user specified positionalconstraints. That is, the user selects a small set of points in thesource domain which is visualized as a textured image, and interactivelydrags them around to new positions.

FIG. 1 shows a schematic flowchart of a method 100 for image deformationaccording to an embodiment of the invention. A user enters a set ofpositional constraints 110. In step 140, an image deformation map isdetermined, based on the positional constraints and a domain, i.e. ashape, extracted from an image 120 in step 130. In step 140, the imagedeformation map is applied to the extracted shape, in order to obtain adeformed image 170. The deformed image 170 is output in step 160.

According to the present embodiment of the invention, the imagedeformation map is a harmonic planar map.

A real-valued function u(x,y): Ω→R, is harmonic if it satisfies theLaplace equation u_(xx)+u_(yy)=0. A planar harmonic map is a map ƒ: Ω→R²where ƒ=[u(x,y),v(x,y)], and the component functions u and v areharmonic. When R² is identified with the complex plane C using complexnotation z=x+iy, one can express f as a complex-valued functionƒ(z)=u(z)+iv(z). Thus a complex valued harmonic function ƒ: Ω⊂C→C can beinterpreted as a harmonic planar map. A complex-valued function ƒ=u+ivis said to be holomorphic if it satisfies the Cauchy-Riemann equations.ƒ(z) is called anti-holomorphic if ƒ(z) is holomorphic.

The Wirtinger derivatives of a complex-valued function are defined asfollows: ƒ _(z) :=½(ƒ_(x)−iƒ_(y)) and ƒ _(z) :=½(ƒ_(x)+iƒ_(y)).

Using the Wirtinger derivatives, the Cauchy-Riemann equations areexpressed concisely as ƒ _(z) ≡0. The Laplacian ƒ_(xx)+ƒ_(yy) isconcisely expressed as 4ƒ _(zz,) hence the condition for harmonicity (ofRef and Imf) is simply ƒ _(zz)≡0. It is then easy to see that the realand imaginary components of a holomorphic function are harmonic. Thesame is true for anti-holomorphic functions. The converse is not truesince ƒ _(zz)≡0 does not imply ƒ _(z) ≡0 nor ƒ_(z)≡0.

The Wirtinger derivatives are useful for defining many of the localgeometric distortion quantities of a C¹ planar map. The local change inarea is: det(J_(ƒ))=|ƒ_(z)|²−|ƒ _(z) |², where J_(f) is the 2×2 Jacobianmatrix of the map. A C1 map ƒ is locally injective and orientationpreserving at a neighborhood of a point z if det(J_(ƒ)(z))>0, oralternatively if |ƒ_(z)(z)|>|ƒ _(z) (z)|.

Isometric and conformal distortion measures are typically formulated interms of the singular values of J_(f), 0≤σ₂≤σ₁:σ₁(z)=|ƒ_(z)|+|ƒ _(z) |,σ₂(z)=∥ƒ_(z)|−|ƒ _(z) ∥  (2)For orientation preserving locally injective maps, the expression on theright simplifies to σ₂(z)=|ƒ_(z)|−|ƒ _(z) |>0. The little dilatationk(z)=|ƒ _(z) (z)|/|ƒ_(z)(z)| is a measure of conformal distortion.

The invention generalizes existing techniques that are designed forsimply-connected domains to multiply-connected domains. The constructionis based on the following theorem, proven by the inventors:

Theorem (Harmonic Decomposition). Let Ω be a multiply connected planardomain with N holes K₁, . . . , K_(N), and choose N arbitrary pointsρ_(i) in K_(i). Then, any harmonic map f: Ω→C can be represented as:

$\begin{matrix}{{f(z)} = {{\overset{\sim}{\Phi}(z)} + \overset{\_}{\overset{\sim}{\Psi}(z)} + {\sum\limits_{i = 1}^{N}{\omega_{i}{Ln}{{z - \rho_{i}}}}}}} & (3)\end{matrix}$where {tilde over (Φ)}, {tilde over (Ψ)}: Ω→C are holomorphic functions,and ω₁, . . . , ω_(N) are some complex coefficients.

The harmonic decomposition (3) is unique up to an additive complexconstant that can be chosen by setting e.g. {tilde over (Ψ)}(z₀)=0 foran arbitrary point z0∈Ω. Theorem Harmonic Decomposition extends asimilar result for simply-connected domains for which the summation termin (3) is missing. While {tilde over (Φ)}(z)+{tilde over (Ψ)}(z) isstill harmonic on multiply-connected domains, the additional summationterm is essential to represent certain harmonic maps and without it, therepresentation is incomplete.

The main theoretical result in [Chen and Weber 2015] is a boundary valuecharacterization of the injectivity and the conformal and isometricdistortion of planar harmonic maps on simply-connected domains. Wegeneralize this result to the multiply-connected case:

Theorem (Bounded Distortion). A planar harmonic map ƒ: Ω→C on amultiply-connected domain with exterior boundary curve γ0 orientedcounterclockwise and interior boundary curves γ₁ . . . γ_(N) orientedclockwise, is locally injective with an upper bound k∈[0, 1) on theconformal distortion, a lower bound σ2>0 on the small singular value ofthe Jacobian, and an upper bound σ1<∞ on the large singular value atevery point z in Ω if and only if:

$\begin{matrix}{{{{\oint_{\gamma_{0}}{\frac{f_{z}^{\prime}(w)}{f_{z}(w)}{dw}}} + {\sum\limits_{i = 1}^{N}{\oint_{\gamma_{i}}{\frac{f_{z}^{\prime}(w)}{f_{z}(w)}{dw}}}}} = 0},} & \left( {4a} \right) \\{{0 \leq {k(w)} \leq {k\mspace{14mu}{\forall{w \in {\partial\Omega}}}}},} & \left( {4b} \right) \\{{{\sigma_{1}(w)} \leq {\sigma_{1}\mspace{14mu}{\forall{w \in {\partial\Omega}}}}},} & \left( {4c} \right) \\{\sigma_{2} \leq {{\sigma_{2}(w)}\mspace{14mu}{\forall{w \in {{\partial\Omega}.}}}}} & \left( {4d} \right)\end{matrix}$Intuitively, Theorem Bounded Distortion states that in order to induceglobal distortion bounds, it is sufficient to bound the distortion of aharmonic map on the boundary curves γ_(i) as long as (4a) holds. Theboundary integral condition (4a) is equivalent to the condition that ƒzis nonvanishing throughout the domain. Namely that ƒ_(z)(z)≠0. Incontrast to the simply-connected domain, the integral over the exteriorboundary γ₀ may be nonzero. This gives rise to maps which are notregular homotopic to the identity map but can still be locallyinjective.

FIG. 2 shows a schematic drawing of a knot. The left side shows adoubly-connected domain Ω. On the right side, an image of a map f: Ω→Cis shown. While f is locally injective, f and the identity are notregular homotopic. Intuitively, one cannot “morph” the annulus on theleft to the one on the right continuously without breaking localinjectivity. The conditions of Theorem 4.2 on f hold. The integral(Equation (4a)) over the exterior boundary is 2π i while the integralover the interior one is −2π i.

As opposed to [Chen and Weber 2015], the inventive deformation algorithmdoes not explicitly impose user-defined bounds k, σ₁, σ₂ on thedistortion. Instead, an overall lower average distortion is favored. Byusing energies that become infinite when the map degenerates at aboundary point, a finite bound on distortion is naturally formed on theboundary. Theorem Bounded Distortion then assures that this bound isalso global. At each iteration, the inventive algorithm certifies themap as locally injective by verifying that (4a) holds, as well as thefollowing:|ƒ_(z)(w)|>|ƒ _(z) (w)| ∀w∈∂Ω.  (5)

How to enforce (4a) and (5) in practice is explained in the followingdescription.

According to the present embodiment, equation (3) is discretized toobtain a finite dimensional space of harmonic maps. The holomorphicfunctions {tilde over (Φ)} and {tilde over (Ψ)} are discretized by usingthe Cauchy complex barycentric coordinates [Weber et al, 2009] which arederived by approximating the boundary of the domain using a polygonalshape. Let {circumflex over (P)}={z₁, z₂, . . . , z_(m)}, z_(j)∈C be thevertices of a multiply-connected planar polygon (the cage). The verticesof the cage are split into sets such that the first set represents theexterior polygon oriented counterclockwise while the following setsrepresent the interior polygons oriented clockwise.

FIG. 3 shows notations. The purple area is a doubly-connected domain Ωbounded by a polygon P with 11 vertices. The solid black polygon is anoutward offset cage {circumflex over (P)}=γ₀∪γ₁ with vertices z_(j)(cyan) composed of exterior polygon γ₀ oriented counterclockwise, andinterior clockwise polygon γ₁. The point z (green) is in Ω or on ∂Ω.A_(j)=z_(j)−z_(j−1), B_(j)(z)=z_(j)−z, and [v_(i), v_(i+1)] is a smallsegment (orange) on P. The point ρ1 (red) is in γ₁.

The holomorphic functions {tilde over (Φ)} and {tilde over (Ψ)} arerepresented as:

$\begin{matrix}{{\overset{\sim}{\Phi}(z)} = {\sum\limits_{j = 1}^{m}{{{\overset{\sim}{C}}_{j}(z)}\varphi_{j}\mspace{14mu}{\overset{\sim}{\Psi}(z)}{\sum\limits_{j = 1}^{m}{{{\overset{\sim}{C}}_{j}(z)}\psi_{j}}}}}} & (6)\end{matrix}$where {tilde over (C)}_(j(z)) is the j^(th) holomorphic Cauchybarycentric coordinate associated with vertex z, and φ_(j), ψ_(j) arecomplex coefficients.

{tilde over (C)}_(j(z)) possess a rather simple closed-form expression(see Weber Ofir. 2010. Hybrid Methods for Interactive ShapeManipulation. Ph.D. Dissertation. Technion—Israel Institute ofTechnology). While Weber et al. [2009] assumed that {circumflex over(P)} is simply connected, their derivation can be extended to the caseof multiple boundary components due to the fact that Cauchy's integralformula still holds. This can be done by integrating the Cauchy kernelover all the boundary components. Practically, representing holomorphicfunctions on such a domain simply boils down to using additional basisfunctions that correspond to vertices that lie on the interiorboundaries. The complex barycentric map induced by the Cauchycoordinates on multiply-connected domains has a nonempty null space of(complex) dimension 2N which corresponds to a similarity transformationof the vertices of the interior boundary polygons. To remove thisambiguity, the first two complex coefficients of each hole are fixed tozero.

The holomorphic functions {tilde over (Φ)} and {tilde over (Ψ)} in theharmonic decomposition (3) are substituted with the expressions from(6).

The additional term Σ_(i=1) ^(N)Ln|z−ρ_(i)| in (3) is substituted withthe expression:

$\begin{matrix}{{\sum\limits_{j = {m + 1}}^{m + N}{\left( {\varphi_{j} + {\overset{\_}{\psi}}_{j}} \right){Ln}{{z - \rho_{j - m}}}}},} & (7)\end{matrix}$where φ_(j), ψ_(j) are 2N additional complex coefficients. Splittingω_(i) into two variables creates some unnecessary redundancy,nevertheless it greatly simplifies the exposition as well as theimplementation.

Optionally, this redundancy may be removed by enforcing a constraintφ_(j)=

, j=m+1, . . . , m+N for each hole. Finally, rearranging terms anddenoting n=m+N leads to:

$\begin{matrix}{{{f(z)} = {{\sum\limits_{j = 1}^{n}{{C_{j}(z)}\varphi_{j}}} + \overset{\_}{\sum\limits_{j = 1}^{n}{{C_{j}(z)}\psi_{j}}}}}{{C_{j}(z)} = \left\{ {\begin{matrix}{\overset{\sim}{C}}_{j{(z)}} & {{j = 1},\ldots\mspace{14mu},m} \\{{Ln}{{z - \rho_{j - m}}}} & {{j = {m + 1}},\ldots\mspace{14mu},n}\end{matrix}.} \right.}} & (8)\end{matrix}$

According to the invention, Equation (8) is the closed-form expressionfor a finite dimensional harmonic subspace with 2n complex variables{φ_(j),

_(j)}_(j=1) ^(n). The dimension of the (complex) null space is 4N+N+1,where 4N is due to the 2 complex DOF per hole for each of the two Cauchybarycentric maps (as explained above), N is due to the redundancyintroduced by Equation (7) and the term 1 is due to the constant DOF ofthe harmonic decomposition (3).

By linearity of the Wirtinger operators, the fact that the derivative ofthe Cauchy coordinates with respect to z is 0, and assuming that φ_(j)=

, ∀j=m+1, . . . , n the Wirtinger derivatives of the harmonic map, justlike the map itself, can be expressed in closed-form:

$\begin{matrix}{{{f_{z}(z)} = {\sum\limits_{j = 1}^{n}{{D_{j}(z)}\varphi_{j}}}},{{\overset{\_}{f_{\overset{\_}{z}}}(z)} = {\sum\limits_{j = 1}^{n}{{D_{j}(z)}\psi_{j}}}},{{D_{j}(z)} = \left\{ \begin{matrix}{{\overset{\sim}{D}}_{j}(z)} & {{j = 1},\ldots\mspace{14mu},m} \\\frac{1}{z - p_{j - m}} & {{j = {m + 1}},\ldots\mspace{14mu},n}\end{matrix} \right.}} & (9)\end{matrix}$

The derivatives of the Cauchy barycentric coordinates {tilde over(D)}_(j)(z) are holomorphic and so is

$\frac{1}{z - p_{j - m}},$hence it is clear that f_(z) and ƒ_(z) are both holomorphic (though notnecessarily integrable).

The main approach in interactive shape deformation is to design andminimize a distortion energy that aggregates a differential quantitythat strives to keep the map as-isometric-as-possible. Such differentialquantities are minimized when the map is locally isometric. However,since a perfect isometry cannot be obtained in general in the presenceof positional constraints, the particular definition of proximity toisometry greatly affects the end result. It is a common practice tomeasure distortion at a point z as a function E(z)=E(σ₁(z), σ₂(z)) ofthe two singular values of the Jacobian matrix J_(f)(z). Such a functionis invariant to rigid motions of both the source and the target domains.E(z) should be minimized (at a single point z) if σ₁(z)=σ₂(z)=1.Energies that minimize conformal distortion are also popular. Ourharmonic subspace contains many pure (with zero distortion) conformalmaps, hence, if desired, we can simply restrict the subspace byeliminating the ψ_(j) variables and use any of our isometric energies toregularize the result.

FIG. 4 shows harmonic deformations of a triply-connected domain. Thedomain on the left has N=2 holes in it. The result in the middleminimizes the Advanced MIPS energy, while the one on the right is a pureconformal map with least isometric distortion.

Appropriate choices of isometric energy for designing locally injectivemaps become infinite when the map collapses locally, which serves as anatural barrier term. Since our optimization depends on high orderderivatives, we focus on energies which are smooth in our variablesφ_(j), ψ_(j). To this end, the distortion measure at z is expressed as asmooth function of |f_(z)|² and |ƒ _(z) |². Since |f_(z)|² and |ƒ _(z)|² are quadratic functions in φ_(j), ψ_(j), the composition is alsosmooth.

FIG. 5 (Table 1) lists several possible choices for such smoothisometric energies. These include the exponential symmetric Dirichletisometric energy: E_(exp)=exp(s·E_(iso)) (Rabinovich Michael, PoranneRoi, Panozzo Daniele, and Sorkine-Hornung Olga. 2017. Scalable LocallyInjective Mappings. ACM Transactions on Graphics 36, 2, Article 16(2017)). In the first column, an isometric measure is shown in terms ofthe singular values. In the second column, the isometric measure isshown in terms of |f_(z)|², |fƒ _(z) |², which are denoted x, y forclarity. It is straightforward to see that these expressions are smoothin x, y wherever x≠y, which is always the case for the locally injectivemaps according to the invention. The third column shows the parametersthat are needed for instantiating the gradient (17) and the Hessian (18)expressions. The forth column shows the four eigenvalues of K (notnecessarily sorted). Potentially negative ones are highlighted.

The motivation to use E_(exp) is to penalize more drastically highvalues of the distortion measure.

FIG. 6 shows a comparison of E_(iso) (left) and E_(exp) (right). Astraight bar shape is deformed into an extreme pose. Exponentiating theE_(iso) distortion measure reduces the maximal distortion (red color)but increases the average distortion. Both color maps show E_(iso). Thisallows the user to trade-off low average versus low maximal distortion.Finally, the Advanced MIPS energy (Fu Xiao-Ming, Liu Yang, and GuoBaining. 2015. Computing locally injective mappings by advanced MIPS.ACM Transactions on Graphics 34, 4 (2015), 71) is included, whichprovides user-controlled balance between area preservation and anglepreservation.

According to the invention, the isometric distortion of a planarharmonic map f is defined as a boundary integral over the pointwisedistortion quantity E(w):E ^(ƒ)=

_(∂Ω) E(w)ds.  (12)

The superscript f is used to denote the boundary integrated distortionand omitted in order to denote the pointwise quantity. Theorem BoundedDistortion ensures that for any bound on E(w) along the boundary, aglobal upper bound on σ₁(z), and a global lower bound on σ₂(z) exist.Hence, a global bound on E(z) is naturally formed. Alternatively, anarea integral may be chosen instead of the boundary one.

The gradient and the Hessian of the overall isometric energy are:∇E ^(ƒ)=

_(∂Ω) ∇E(w)ds.  (13)∇² E ^(ƒ)=

_(∂Ω)∇² E(w)ds.  (14)

The gradient and the Hessian of our energy measures have relativelysimple closed-form expressions. Let D=D1, D2, . . . , Dn)∈C^(1×n) be acomplex row vector, where D_(j) is defined as in (9). Bold symbols areused to denote real vectors and matrices. Define the real matrixD(notethe bold symbol) as:

$\begin{matrix}{D = {\begin{bmatrix}{{Re}(D)} & {- {{Im}(D)}} \\{{Im}(D)} & {{Re}(D)}\end{bmatrix} \in \Re^{2{x2n}}}} & (15)\end{matrix}$

The complex Wirtinger derivatives are expressed as 2×1 real vectors:

$\begin{matrix}{{f_{z} = \begin{bmatrix}{{Re}\left( f_{z} \right)} \\{{Im}\left( f_{z} \right)}\end{bmatrix}},{\overset{\_}{f_{z}} = \begin{bmatrix}{{Re}\left( \overset{\_}{f_{\overset{\_}{z}}} \right)} \\{{Im}\left( \overset{\_}{f_{\overset{\_}{z}}} \right)}\end{bmatrix}}} & (16)\end{matrix}$

The gradient of E with respect to the 4n real variables is:

$\begin{matrix}{{\nabla{E(z)}} = {{2\begin{bmatrix}{\alpha_{1}D^{T}f_{z}} \\{\alpha_{2}D^{T}\overset{\_}{f_{\overset{\_}{z}}}}\end{bmatrix}} \in \Re^{4{nx}\; 1}}} & (17)\end{matrix}$

where α₁, α₂ are real parameters which depend on the particular choiceof energy. Table 1 provides the parameters needed to instantiate someparticular energies (out of many possible). The expression for the 4n×4nHessian of E(z) at a single point is:

$\begin{matrix}{{\nabla^{2}{E(z)}} = {{\begin{bmatrix}D^{T} & 0 \\0 & D^{T}\end{bmatrix}{K\begin{bmatrix}D & 0 \\0 & D\end{bmatrix}}} \in \Re^{4{nx}\; 4n}}} & (18)\end{matrix}$where K is the a 4×4 real matrix:

$\begin{matrix}{K = {\begin{bmatrix}{{2\alpha_{1}I} + {4\beta_{1}f_{z}f_{z}^{T}}} & {4\beta_{3}f_{z}{\overset{\_}{f_{\overset{\_}{z}}}}^{T}} \\{4\beta_{3}\overset{\_}{f_{\overset{\_}{z}}}f_{z}^{T}} & {{2\alpha_{2}I} + {4\beta_{2}\overset{\_}{f_{\overset{\_}{z}}}{\overset{\_}{f_{\overset{\_}{z}}}}^{T}}}\end{bmatrix} \in \Re^{4x\; 4}}} & (19)\end{matrix}$

The parameters β₁, β₂, β₃ are also energy dependent (see Table 1).

The point-to-point (P2P) metaphor is an intuitive drag-and-dropuser-interface that best fits interactive deformation tasks. The user ofthe inventive method can add or remove (by clicking) positionalconstraints at any time during interaction. Dragging the P2P handlessignals the application to invoke the optimization and to render theupdated result. The P2P constraints are incorporated as softconstraints. The P2P energy is defined as:

${E_{p\; 2p}^{f} = {\frac{1}{2}{\sum\limits_{i = 1}^{P}\;{{{f\left( p_{i} \right)} - q_{i}}}^{2}}}},$where p_(i)∈P⊂Ω is the position of the i^(th) handle in the sourcedomain, and q^(i)∈C is its target desired position. The gradient and theHessian of E_(p2p) ^(ƒ) are derived as follows.

The P2P energy at a single handle can be expressed using complexnumbers:E _(p2p)(p _(i))=½|ƒ(p _(i))−q _(i)|²=½|C(p _(i))φ+

−q _(i)|²,  (F.1)where C(p_(i))=(C₁(pi), C₂(pi), . . . , C_(n)(pi))∈C^(1×n) is a complexrow vector, where is defined in (8), and φ,ψ∈C^(n×1). The completeenergy is given by the sum:

$\begin{matrix}{E_{p\; 2p}^{f} = {\sum\limits_{i = 1}^{P}\;{{E_{p\; 2p}\left( p_{i} \right)}.}}} & \left( {F{.2}} \right)\end{matrix}$

Equation (F.1) can be expressed using real numbers:

${{E_{p\; 2p}\left( p_{i} \right)} = {\frac{1}{2}{{{{Q\left( p_{i} \right)}\begin{bmatrix}\varphi \\

\end{bmatrix}} - \begin{bmatrix}{{Re}\left( q_{i} \right)} \\{{Im}\left( q_{i} \right)}\end{bmatrix}}}^{2}}},$where Q(p_(i)) is a real matrix defined as:

[ Re ( C ( p i ) ) - Im ( C ( p i ) ) Re ( C ( p i ) ) - Im ( C ( p i )) Im ( C ( p i ) ) Re ( C ( p i ) ) - Im ( C ( p i ) ) - Re ( C ( p i )) ] ∈ 2 × 4 ⁢ n

E_(p2p)(p_(i)) is a quadratic form with constant PSD Hessian matrix:∇₂ E _(p2p)(p _(i))=Q ^(T)(p _(i))Q(p _(i))∈

^(4n×4n)  (F.3)

The gradient of the quadratic form is:

$\begin{matrix}{{\nabla^{2}{E_{p\; 2p}\left( p_{i} \right)}} = {{{Q^{T}\left( p_{i} \right)}{{Q\left( p_{i} \right)}\begin{bmatrix}\varphi \\

\end{bmatrix}}} - {{{Q^{T}\left( p_{i} \right)}\begin{bmatrix}{{Re}\left( q_{i} \right)} \\{{Im}\left( q_{i} \right)}\end{bmatrix}}.}}} & \left( {F{.4}} \right)\end{matrix}$

In order to sum the gradients and the Hessians over all the P2P handlesP, one constructs a matrix Q by stacking Q(p_(i)):

Q = [ Q ⁡ ( p 1 ) Q ⁡ ( p 2 ) ⋮ Q ( p  P  ] ∈ 2 ⁢  P  × 4 ⁢ n ( F ⁢ .5 )

Finally, the full Hessian and full gradient are given by:

$\begin{matrix}{{\nabla^{2}E_{p\; 2p}^{f}} = {Q^{T}Q}} & \left( {F{.6}} \right) \\{{\nabla E_{p\; 2p}^{f}} = {Q^{T}\left( {{Q\begin{bmatrix}\varphi \\

\end{bmatrix}} - q} \right)}} & \left( {F{.7}} \right)\end{matrix}$where q is a real column vector stacking the target positions:q=└Re(q ₁),Im(q ₂), . . . ,Re(q _(|P|)),Im(q _(|P|))┘^(T)∈

^(2|P|×1)  (F.8)

The full energy of the unconstrained minimization problem is:E _(Def) ^(ƒ) =E ^(ƒ) +λE _(p2p) ^(ƒ),  (20)where λ is a user-defined weight that balances the two terms.

The Hessian of E_(p2p) ^(ƒ) is positive semi-definite (PSD). However,the Hessian of E^(f) is, in general, not PSD, and neither is the Hessianof E_(Def) ^(ƒ) (Equation (20)). The key observation here is that theclosest (in Frobenius norm) PSD matrix to the Hessian ∇²E(z)∈

^(4n×4n) at a single point z, can be expressed in closed-form. With thatat hand Equation (14) is substituted with:∇² E ^(ƒ+)=

_(∂Ω)∇² E ⁺(w)ds,  (21)where ∇₂E⁺(w) is the closest PSD matrix to ∇²E(w). Since, the integral(or sum) of PSD matrices is PSD (due to the convex cone structure of thePSD set), it is clear that the modified Hessian of the isometric energy(21) is guaranteed to be PSD. Further, the dimension of the null spaceof (21) in the case of E_(iso) and E_(exp) is 2 or 3. Therefore, with 2or more positional constraints, ∇²E^(ƒ+) is nonsingular.

This application refers to the present variant of Newton's method thatcomputes the closest PSD matrix to ∇²E^(ƒ) as Newton-Eigen. The modifiedNewton iterations are dramatically faster to compute. Moreover,throughout extensive experiments, the inventors observed that the localmodification leads to iterations which are more effective, where lessiterations are required for convergence (FIGS. 9, 12).

Computing the eigen-factorization of the local Hessian matrix ∇²Enumerically is impractical. In order to derive an analyticeigen-factorization of a 4n×4n matrix, it is first observed that each ofthe nontrivial eigenvectors (with nonzero eigenvalue) of ∇²E can beexpressed as a product of

$B = \begin{bmatrix}D & 0 \\0 & D\end{bmatrix}$and a corresponding eigenvector of K (Equation (19)). This is due to thefact that B has 4 rows that are orthogonal to each other (easy to verifyby computing the inner product of each pair) and all the rows have thesame norm. Furthermore, the eigenvalues of K and ∇2E are the same up toa positive scale. Hence, the above problem can be reduced toanalytically computing the eigenvalues of a 4×4 matrix. This is stillchallenging as these are the roots of a 4th order polynomial. Denote the(unsorted) eigenvalues of K as (λ1, λ2, λ3, λ4). It is further observedthat λ₁=2α1 and λ₂=2α2. To see why 2α₁ is an eigenvalue, subtract 2α₁from the diagonal of K. The first two rows of K−2α₁I are:ƒ_(z)└4β₁ƒ_(z) ^(T) 4β₃ ƒ _(z) ^(T),  (22)where we can see that the 1×4 row vector on the right hand side of (22)is multiplied by two scalars (the elements of f_(z)), hence, these tworows are linearly dependent, meaning that the matrix K−2α₁I is singularand 2α₁ is a root of the characteristic polynomial. Similarly, it can beshown that 2α₂ is an eigenvalue. Knowing that 2α₁ and 2α₂ are twoeigenvalues of K, the other two eigenvalues one can compute by directlysolving the quartic characteristic equation. The derivation is long butstraightforward, hence omitted. One obtains these expressions:λ_(3,4) =s ₁±√{square root over (s ₂ ²+16β₃ ²|ƒ_(z)|²|ƒ _(z) |²)}.  (23)where s_(1,2)=α₁+2β₁|ƒ_(z)|²±(α₂+2β₂|ƒ _(z) |²).

The corresponding four eigenvectors are:

$\begin{matrix}{\left( {{{Im}\left( f_{z} \right)},{- {{Re}\left( f_{z} \right)}},0,0} \right)\left( {0,0,{{Im}\left( \overset{\_}{f_{\overset{\_}{z}}} \right)},{- {{Re}\left( \overset{\_}{f_{\overset{\_}{z}}} \right)}}} \right)\left( {{{Re}\left( f_{z} \right)},{{Im}\left( f_{z} \right)},{t_{1}{{Re}\left( f_{\overset{\_}{z}} \right)}},{t_{1}{{Im}\left( f_{\overset{\_}{z}} \right)}}} \right)\left( {{{Re}\left( f_{z} \right)},{{Im}\left( f_{z} \right)},{t_{2}{{Re}\left( f_{\overset{\_}{z}} \right)}},{t_{2}{{Im}\left( f_{\overset{\_}{z}} \right)}}} \right){{where}\text{:}}{t_{1,2} = {\frac{\lambda_{3,4} - {2\;\alpha_{1}} - {4\;\beta_{1}{f_{\overset{\_}{z}}}^{2}}}{4\;\beta_{3}{f_{\overset{\_}{z}}}^{2}}.}}} & (24)\end{matrix}$

The signs of these eigenvalues depend on the particular choice ofisometric energy, therefore in the most general case, the 4 eigenvaluesare evaluated and negative ones are substituted with 0 to obtain themodified Hessian. For particular energy choices, it is possible tosimplify matter even more. For the first two energies listed in Table 1,only λ₁=2α₁ can be (at times) negative. This allows to directly expressthe modified matrix. For example, for E_(iso), it turns out that λ₃, λ₄have quite simple expressions:λ₃=4(1+3(|ƒ_(z)|+|ƒ _(z) |⁻⁴),λ₄=4(1+3(|ƒ_(z)|−|ƒ _(z) |⁻⁴).

The spectrum is then sorted as follows: 2α₁≤λ₃≤2α₂≤λ₄. With thisinformation at hand, the modification is done by checking for the signof α₁, and if it is negative, K is substituted in (18) with:

$\begin{matrix}{K^{+} = \begin{bmatrix}{\left( {\frac{2\;\alpha_{1}}{{f_{z}}^{2}} + {4\;\beta_{1}}} \right)f_{z}f_{z}^{T}} & {4\;\beta_{3}f_{z}{\overset{\_}{f}}_{\overset{\_}{z}}^{T}} \\{4\;\beta_{3}{\overset{\_}{f}}_{\overset{\_}{z}}^{T}f_{z}^{T}} & {{2\;\alpha_{2}I} + {4\;\beta_{2}\overset{\_}{f_{\overset{\_}{z}}}{\overset{\_}{f}}_{\overset{\_}{z}}^{T}}}\end{bmatrix}} & (25)\end{matrix}$

To ensure that the map is locally injective at each iteration, it needsto be verified that Conditions (4a) and (5) hold, otherwise one needs tobacktrack during line search. As (5) involves infinite number ofinequalities: |ƒ_(z)(w)|>|ƒ _(z) (w)| ∀w∈∂Ω, the method uses a simplersufficient condition based on a finite number of conditions, utilizingthe fact that the Wirtinger derivatives are Lipschitz continuous.

Condition (5) can be enforced on the entire boundary by enforcing it(individually) on many small boundary segments. For each segment[v_(i),v_(i+1)] (see FIG. 4 for notations), one computes lower and upperbounds for |ƒ_(z)| and |ƒ _(z) | respectively as follows:

${{f_{z}}_{\min}:={{\frac{{{f_{z}\left( v_{i} \right)}} + {{f_{z}\left( v_{i + 1} \right)}}}{2} - \frac{L_{f_{z}}l}{2}} \leq {\min\limits_{w \in {\lbrack{v_{i},v_{i + 1}}\rbrack}}{{f_{z}(w)}}}}},{{f_{\overset{\_}{z}}}_{\max}:={{\frac{{{f_{\overset{\_}{z}}\left( v_{i} \right)}} + {{f_{\overset{\_}{z}}\left( v_{i + 1} \right)}}}{2} + \frac{L_{f_{\overset{\_}{z}}}l}{2}} \leq {\max\limits_{w \in {\lbrack{v_{i},v_{i + 1}}\rbrack}}{{{f_{\overset{\_}{z}}(w)}}.}}}}$

L_(ƒ) _(z) and L_(ƒ) _(z) are the corresponding Lipschitz constants onthe segment, and I=|v_(i)−v_(i+1)|. Condition (5) is then substitutedwith:|ƒ_(z)(v _(i))|−|ƒ _(z) (v _(i))|+|ƒ_(z)(v _(i+1))|−|ƒ _(z) (v_(i+1))|≥(L _(ƒ) _(z) +L _(ƒ) _(z) )l,  (26)or more concisely:σ₂(v _(i))+σ₂(v _(i+1))≥(L _(ƒ) _(z) +L _(ƒ) _(z) )l.  (27)

It is crucial that the sufficient condition above is as tight aspossible, as otherwise the Newton step size may become too small andwill stop the Newton iterations prematurely. The condition becomestighter as the Lipschitz constants and/or the length of the segment 1decrease. Since denser sampling requires more computations, it isadvised to obtain as small as possible Lipschitz constants.

According to the invention, the following Lipschitz constants may beused:

$L_{f_{z}} = {\frac{{{f_{z}^{\prime}\left( v_{i} \right)}} + {{f_{z}^{\prime}\left( v_{i + 1} \right)}}}{2} + {\frac{l}{2}\left( {{\sum\limits_{j = 1}^{m}\;{L_{j}{{s_{j} - s_{j + 1}}}}} + {\sum\limits_{j = {m + 1}}^{n}\;{L_{j}^{h}{\varphi_{j}}}}} \right)}}$L f z _ =  f z _ ′ ⁡ ( v i )  +  f z _ ′ ⁡ ( v i + 1 )  2 + l 2 ⁢ ( ∑ j= 1 m ⁢ ⁢ L j ⁢  t j - t j + 1  + ∑ j = m + 1 n ⁢ ⁢ L j h ⁢  j  )where

${L_{j} = \frac{1}{2\pi\;{d^{\; 2}\left( z_{j} \right)}}},{L_{j}^{h} = \frac{2}{d^{3}\left( {\rho_{j} - m} \right)}},$d(z) is the distance from a point z to the segment,

s j = φ j - φ j - 1 z j - z j - 1 ⁢ ⁢ and ⁢ ⁢ t j = j - j - 1 z j - z j - 1.

FIG. 7 shows a comparison of the Lipschitz constants obtained with theabove formulas against known methods. More particularly, FIG. 7 shows onthe left a comparison of the Lipschitz constants of fz on 529 boundarysegments. L_(BDHM) is the approach taken in [Chen and Weber 2015], whileL_(ours) is the improved constants. The above constants are smaller onall the segments with an average (number in parentheses above) which is×6.74 smaller. The right side shows a comparison of convergence behaviorof the inventive solver using BDHM constants (black) vs. invention(red). See how the iterations of the Raptor and Horse deformations stopimmaturely when L_(BDHM) is used. In summary, the above given constantsare significantly smaller and lead to faster convergence and increasedrobustness.

Turning now to Condition (4a) requiring that the boundary integral (orequivalently the number of zeros of f_(z)) is 0. The following conditionis applicable to multiply-connected domains, and is both necessary andsufficient as long as (27) holds:

Theorem 8.1. Let g(z) be an L-Lipschitz continuous holomorphic functionin a neighborhood of a simple open curve with length l and two endpointsv_(i),v_(i+1)∈C such that:|g(v _(i))+|g(v _(i+1))|>Ll,  (28)then:

$\begin{matrix}{{\int_{v_{i}}^{v_{i + 1}}{\frac{g^{\prime}(z)}{g(z)}{dz}}} = {{\ln{\frac{g\left( v_{i + 1} \right)}{g\left( v_{i} \right)}}} + {i\;{Arg}{\frac{g\left( v_{i + 1} \right)}{g\left( v_{i} \right)}}}}} & (29)\end{matrix}$where Arg is the principle branch of the complex argument function.

Theorem 8.1 is applicable to f_(z) along any line segment[v_(i),v_(i+1)] since the inventive algorithm always verifies that (27)holds (which implies (28)). Consequently, under the assumption that (27)holds, f_(z) does not vanish inside the multiply-connected polygon P ifand only if:

$\begin{matrix}{{{\sum\limits_{j = 0}^{N}{\sum\limits_{i}^{\;}{{Arg}\frac{f_{z}\left( v_{i + 1}^{j} \right)}{f_{z}\left( v_{i}^{j} \right)}}}} = 0},} & (30)\end{matrix}$where the first summation is over the polygonal loops, and the secondone is over the segments in each loop.

To conclude, the method first verifies that (27) holds on all theboundary segments and if so, the sum in (30) is simply computed. If itis zero, the map is guaranteed to be locally injective everywhere.

In general, Newton's method searches for a stationary point of a scalarfunction E(x) of a vector variable x∈R^(n) by iteratively approximatingit using its second order Taylor expansion:E(x)=E(x _(k) +Δx)≈E(x _(k))+Δx ^(T) ∇E(x _(k))+½Δx ^(T)∇² E(x _(k))Δxwhere Δx=x−x_(k)∈R^(n) is the search direction, and ∇E∈R^(n) and∇²E∈R^(n×n) denotes the gradient and the Hessian of E respectively. Thissecond order polynomial in Δx has vanishing derivative when thefollowing linear system is satisfied:∇² E(x _(k))Δx=−∇E(x _(k))  (1)

The direction Δx decreases the energy only if the second orderpolynomial is convex, or equivalently, the Hessian matrix ∇²E ispositive definite. Hence, when needed, the Hessian is modified. Thesolution to the linear system produces the update for the currentiteration x_(k+1)=x_(k)+Δx. To ensure convergence (to a local minimum)from an arbitrary starting point, a line search along the searchdirection Δx should be used to find a step size t that satisfies theWolfe conditions.

FIG. 8 shows, in pseudo-code, an algorithm (Algorithm 1) for a Newtonsolver used to determine an image deformation map according to anembodiment of the invention.

While the map and its differentials possess closed-form expressions, thecontour integrals that arise in Equations (12), (13), (14) may not besolved analytically; therefore, the method resorts to numericalintegration. Let G={w₁, w₂, . . . , w_(|G|)}⊂∂Ω an be a set of uniformlydistributed samples. Applying the trapezoidal rule to (12), it boilsdown to a simple sum, due to the uniform sampling:

$\begin{matrix}{E^{f} \approx {\frac{1}{G}{\sum\limits_{i = 1}^{G}{E\left( w_{i} \right)}}}} & (31)\end{matrix}$

The gradient of E^(f) is approximated in the same fashion:

$\begin{matrix}{{{\nabla E^{f}} \approx {\frac{1}{G}{\sum\limits_{i = 1}^{G}{\nabla{E\left( w_{i} \right)}}}}},} & (32)\end{matrix}$

However, it turns out that for the integral approximation of theHessian, much coarser approximation can be used with negligibleinfluence on the number of Newton iterations (see FIG. 11). Hence, themethod uses an additional set H of samples which is a subset of G butcontains significantly less samples. The Hessian is approximated:

$\begin{matrix}{{{\nabla^{2}E^{f}} \approx {\frac{1}{H}{\sum\limits_{i = 1}^{H}{\nabla^{2}{E^{+}\left( w_{i} \right)}}}}},} & (33)\end{matrix}$

where ∇²E⁺(w_(i)) is the analytic modification of the pointwise Hessian.We note on a theoretical level, that for any number |H| (as small as itmay be), the Newton step produces a descent direction since the matrix∇²E^(f) in (33) is guaranteed to be PSD (strict positive definitenesscan be ensured by adding E to the diagonal). In the extreme case that His an empty set, the algorithm degrades to gradient descent. Theinventors observed that in practice, choosing |H|=0.1 |G| does notdegrade the effectiveness of the Newton iterations and provides a goodbalance between the time it takes to perform the numerical integrationand the other steps of the algorithm. This setting is used in all theresults. Moreover, since in practice |H| is much larger than thedimension of our subspace, it is not necessary to add E to the diagonal.

To ensure convergence to a local minimum from an arbitrary startingpoint the method uses simple backtracking [Nocedal Jorge and WrightStephen. 2006. Numerical Optimization. Springer New York, Algorithm 3.1]to find a step size that sufficiently decreases the energy (Algorithm 1,Line 17) and in addition ensures that the new map will be locallyinjective using the above-described strategy (Algorithm 1, Line 18).This requires evaluations of the energy and the Wirtinger derivatives atdifferent step sizes.

The algorithm shown in FIG. 8 can be dramatically accelerated byexploiting parallelism. Most of the steps require only basic denselinear algebra operations and are straightforward to implement. Unlikesparse linear algebra operations, the inventive method greatly benefitswhen these operations are performed on the GPU, as will be described inthe following.

The GPU implementation proposed by the invention is designed in a waythat utilizes standard libraries as much as possible. The implementationmostly involves only dense linear algebra operations, while the partsthat require coding of a specialized kernel are relatively simple,embarrassingly parallel operations. The inventors have utilized thecuBLAS (cuBLAS. 2017. CUDA Basic Linear Algebra Subroutines.http://developer.nvidia.com/cublas. (2017)), and cuSolverDN (cuSOLVER.2017. CUDA Linear Solvers Library. http://developer.nvidia.com/cusolver.(2017)) (BLAS and LAPACK implementations) dense linear algebra librariesthat are freely distributed with NVIDIA's CUDA toolkit with doubleprecision accuracy. This not only simplifies the implementation, butalso ensures that the code will be optimized for different GPUarchitectures and will run smoothly on any NVIDIA device. Moreover, asgraphics hardware constantly evolves and these libraries get optimizedfurther, the speed of such implementation is expected to automaticallyimprove over time.

Regarding matrix and vector notations, the Wirtinger derivatives at asingle point can be expressed concisely in vector form:ƒ_(z)(z)=D(z)_(φ),ƒ_(z) =D(z)

,  (G.1)where D(z)=(D₁(z), D₂(z), . . . , D_(n)(z))∈C^(1×n) is a complex rowvector, while φ=(φ₁, φ₂, . . . , φ_(n))^(T) and ψ=(ψ₁, ψ₂, . . . ,ψ_(n))^(T) are complex column vectors.

The complex-valued function f and its differentials have to be evaluatedat a large collection of points inside the domain. It would beconvenient to express these in matrix notations. For example, assumethat one wants to evaluate the map f on a set of points V∈Ω for the sakeof visualization.

In the complex matrix C(V)∈C^(|V|×n) the i^(th) row of C corresponds tothe point v_(i)∈V such that C_(i,j)=C_(j) (v_(i)). Similarly, thecomplex matrix D(V)∈C^(|V|×n) may be used such that D_(i,j)=D_(j)(v_(i)). Using these notations, one may evaluate the map fat the pointsof the set V by computing the matrix product:

$\begin{matrix}{{\underset{\underset{{V}x\; 2}{︸}}{\begin{bmatrix}\Phi & \Psi\end{bmatrix}} = {{C(V)}\underset{\underset{{nx}\; 2}{︸}}{\begin{bmatrix}\varphi & \psi\end{bmatrix}}}},} & \left( {G{.2}} \right)\end{matrix}$followed by the computation of the vector sum Φ+Ψ. Similarly, one canevaluate the Wirtinger derivatives on the set V with the product:

$\begin{matrix}{{\underset{\underset{{V}x\; 2}{︸}}{\begin{bmatrix}f_{z} & \overset{\_}{f_{\overset{\_}{z}}}\end{bmatrix}} = {{D(V)}\underset{\underset{{nx}\; 2}{︸}}{\begin{bmatrix}\varphi & \psi\end{bmatrix}}}},} & \left( {G{.3}} \right)\end{matrix}$

In the preprocessing step (Algorithm 1, line 1), the polygon P istriangulated with a dense triangulation whose vertices form the set V.The triangulation may be used for visualization purposes. Hardwaretexture mapping may be used to render the final deformed image. Therendering can be performed efficiently on the GPU with extremely highmesh- and texture resolutions without affecting the runtime of thesolver. The following matrices are computed and stored on the device'smemory:C(V)∈C ^(|V|×n) C(P)∈C ^(|P|×n)D(G)∈C ^(|G|×n) D(H)∈C ^(|H|×n)  (G.4)where V contains the vertices of the triangulation used for texturemapping, P is the set of positional constraints, G is a dense set ofsamples lying on the polygon P which is used for energy and gradientevaluation, and H is a subset of G used to evaluate the Hessian of theenergy. Once the complex matrix D(H) is constructed, it is converted toreal form D∈R^(2|H|×2n) (where the bold symbol denotes a real matrix) asfollows:

$\begin{matrix}{{D = {\begin{bmatrix}{{Re}\left( D_{1} \right)} & {- {{Im}\left( D_{1} \right)}} \\{{Im}\left( D_{1} \right)} & {{Re}\left( D_{1} \right)} \\{{Re}\left( D_{2} \right)} & {- {{Im}\left( D_{2} \right)}} \\{{Im}\left( D_{2} \right)} & {{Re}\left( D_{2} \right)} \\\vdots & \vdots \\{{Re}\left( D_{H} \right)} & {- {{Im}\left( D_{H} \right)}} \\{{Im}\left( D_{H} \right)} & {{Re}\left( D_{H} \right)}\end{bmatrix} \in R^{2{H}{x2n}}}},} & \left( {G{.5}} \right)\end{matrix}$where the subscripts indicates row numbers such that D_(i) is the i^(th)row of the complex matrix D. Another closely related matrix with thesame size which we denote as {circumflex over (D)} is also stored and isobtained from D by swapping each two pair of rows. That is, the 1^(st)row is swapped with the 2^(nd) one, the 3^(rd) with the 4^(th) and soon.

Next the implementation computes in parallel a matrix with thecoefficients

$L_{i,j} = {{\frac{1}{2\pi\;{d_{i}^{2}\left( z_{j} \right)}}\mspace{14mu}{and}\mspace{14mu} L_{i,j}^{h}} = {\frac{2}{d_{i}^{3}\left( {p_{j} - m} \right)}.}}$These will be used in runtime to compute the Lipschitz constants of theWirtinger derivatives. Unlike the Lipschitz constants, these depend onlyon the cage {circumflex over (P)} and the set G, hence, can beprecomputed. A simple formula for the distance d_(i)(p) from a point pto the i^(th) segment is given in Appendix C of [Chen and Weber 2015].

The last step of the preprocessing computes the matrix Q (see Equation(F.5)) and ∇²E_(p2p) ^(ƒ)=Q^(T)Q, which is a constant matrix sinceE_(p2p) ^(ƒ) is a quadratic form, hence can be computed and stored once.

Next, the computations are described that are needed to be carried outon each Newton iteration. The evaluation of the map f, its Wirtingerderivatives, energy, and gradient are computed using a matrix-matrixmultiplication, where the right term in the product is a “thin” matrix(having many rows but very small number of columns). Such operation isknown to be bandwidth limited, meaning that the computation time isdominated by the memory bandwidth of the GPU device. Hence, whenperforming the product one should strive to minimize the amount ofaccess to the global memory of the device. To this end, the computationare performed using complex numbers and the result is converted back toreal numbers when necessary. During runtime, the Wirtinger derivativesare first evaluated at the samples of G using the product:

$\begin{matrix}{{\underset{\underset{{G}x\; 2}{︸}}{\begin{bmatrix}f_{z} & \overset{\_}{f_{\overset{\_}{z}}}\end{bmatrix}} = {{D(G)}\underset{\underset{{nx}\; 2}{︸}}{\begin{bmatrix}\varphi & \psi\end{bmatrix}}}},} & \left( {G{.6}} \right)\end{matrix}$

Next, the following is being computed in a single dedicated kernel thataccepts the matrix └ƒ_(z) ƒ_(z) ┘ as input. Each thread handles a row of└ƒ_(z) ƒ_(z) ┘ (which contains only two complex numbers). The two realscalars |f_(z)|² and |

|², are computed and used to compute the real scalars α₁, α₂ using(E.3). Finally, the kernel multiplies α₁, α₂ element-wise by thecorresponding matrix entries, and store the result as a matrix [α₁ƒ_(z)α₂

]∈R^(|G|×2), These are embarrassingly parallel operations, which arestraightforward to implement in a kernel, where no synchronization isneeded and the memory access pattern is perfectly regular.

The next step computes the following complex matrix-matrix product:

$\begin{matrix}{{{\underset{\underset{{nx}\; 2}{︸}}{\begin{bmatrix}g_{\varphi} & g_{\psi}\end{bmatrix}} = {\frac{2}{G}{D^{H}(G)}\underset{\underset{{G}x\; 2}{︸}}{\begin{bmatrix}{\alpha_{1}f_{z}} & {\alpha_{2}\overset{\_}{f_{\overset{\_}{z}}}}\end{bmatrix}}}},}\;} & \left( {G{.7}} \right)\end{matrix}$where (⋅)H is the conjugate transpose (hermit) operator. Themultiplication in (G.7) is computed by invoking the gemm cuBLASprocedure. Based on the summation in (32), the (real-valued) expressionof the gradient of the energy is given by:

$\begin{matrix}{{\nabla E^{f}} = {\begin{bmatrix}{Re}_{(G_{\varphi})} \\{Im}_{(G_{\varphi})} \\{Re}_{(G_{\psi})} \\{- {Im}_{(G_{\psi})}}\end{bmatrix} \in {R^{4{nx}\; 1}.}}} & \left( {G{.8}} \right)\end{matrix}$

The gradient of the P2P energy also needs to be computed:

$\begin{matrix}{{\nabla E_{p\; 2p}^{f}} = {{Q^{T}\left( {{Q\begin{bmatrix}\varphi \\\psi\end{bmatrix}} - q} \right)}.}} & \left( {G{.9}} \right)\end{matrix}$

Q∈R^(2|P|×4n) was already computed in preprocessing, and the vectorq∈R^(2|P|×1) contains the target P2P positions (Equation (F.8)).

The computation of the Hessian is the most involved part of theimplementation. First, the 4×4 matrix K_(i) ⁺ that corresponds to thei^(th) sample in H is computed in a dedicated kernel as described above.The real scalars α₁, α₂, β₁, β₂, β₃ (Table 1) are computed. Then, basedon the energy type, the implementation checks for negative eigenvaluesand decides whether to compute K_(i) (Equation (19)) or its modifiedversion K_(i) ⁺. In the case of E_(iso), the matrix K_(i) ⁺ can becomputed directly using (25). Otherwise, equations (23) and (24) areused. The entries of all these 4×4 matrices are stored in a matrixK∈R^(16×|H|) where in each column the 16 values of the correspondingmatrix K_(i) ⁺ are stacked.

Unlike the computation of the gradient, which is bandwidth limited, theassembly of the Hessian is a ALU-bound operation, meaning that itsruntime is governed by the amount of processing units (and their clockspeed) rather than the speed of the global memory. This means that onecan expect greater speedup compared to a CPU implementation. Yet, theamount of arithmetic operations involved is vast. The inventive strategyto use coarser sampling (Equation (33)) greatly reduces the Hessianassembly time.

In order to efficiently implement equation (33), the assembly of theHessian is done separately on each of its four 2n×2n blocks. How tocompute the ∇_(φ) ₂ ²E^(ƒ) block is explained below. The blocks for∇_(φ) ₂ ²E^(ƒ) and ∇_(ψψ) ²E^(ƒ) can be computed analogously. Thesethree blocks are copied into the corresponding position in the 4n×4nfull Hessian matrix. Due to symmetry, there is no need to explicitlycompute ∇_(ψψ) ²E^(ƒ) as it is given by ∇_(ψψ) ²E^(ƒ) ^(T) . Let

$\begin{pmatrix}a_{i} & b_{i} \\c_{i} & d_{i}\end{pmatrix}\quad$be the 2×2 upper-left block of the matrix K_(i) ⁺. Note that K (computedearlier) contains the corresponding values for all samples in its rows.The entries of ∇_(φ) ₂ ²E^(ƒ) can be computed using the followingformula:

$\begin{matrix}{{{\nabla_{\varphi^{2}}^{2}E^{f}} = {\frac{1}{H}{D^{T}\left\lbrack {{RD} + {B\hat{D}}} \right\rbrack}}},} & \left( {G{.10}} \right)\end{matrix}$where R∈R^(2|H|×2|H|) is a diagonal matrix with diagonal entries (a₁,d₁, a₂, d₂, . . . , a_(|H|), d_(|H|)) and B∈R^(2|H|×2|H|) is a diagonalmatrix with diagonal entries (b₁, c₁, b₂, c₂, . . . , b_(|H|), c_(|H|)).The advantage of performing the computation in such a way is the abilityto use the highly-optimized cuBLAS library [cuBLAS 2017]. The productsRD and B{circumflex over (D)} in (G.10) are computed using a diagonalmatrix-matrix multiplication procedure dgmm followed by a matrixsummation (geam) which result in a matrix M=RD+B{circumflex over (D)}.Finally, the result is obtained through a matrix-matrix product D^(T)Musing gemm. This product dominates the runtime of (G.10).

Then, the gradient and the Hessian of E_(p2p) ^(ƒ) are added to those ofE^(f) (with the appropriate user-defined weight λ) to form the gradientand the Hessian of the full energy E^(f) (Equation (20)). The 10N+2redundant real variables that correspond to the dimension of the nullspace of the harmonic subspace are then eliminated as explained above.This is done by removing the corresponding rows and columns from thegradient and the Hessian. After that, the Newton step is computed. Thefreely available cuSolverDN library is used to solve the linear system(1) using the dense GPU-accelerated direct Cholesky factorization.

The solution of the linear system provides the search direction for theNewton iteration. However, this step is not taken immediately. To ensureconvergence to a local minimum which corresponds to a locally injectivemap, a backtracking line search is performed. This requires multipleevaluations of the energy and the local injectivity validation atvarious step sizes, starting from t=1 and dividing it by 2 at each stepuntil the resulted map is accepted. For a map to be accepted, it isverified first that the energy decreases sufficiently (Algorithm 1, Line17) and only then the local injectivity validation (Algorithm 1, Line18) is performed. Both steps require evaluation of the Wirtingerderivatives ƒ_(z) ^(t) and

at all samples, which varies depending on t. In general, evaluatingthese boils down to the same matrix product used in (G.6) for thecomputation of the gradient. Evaluating it many times can easily becomethe algorithm's bottleneck. Luckily, since the derivatives are linear inthe variables, this can achieved with a single matrix product. Thederivatives for step t are given by:

$\begin{matrix}\begin{matrix}{\left\lfloor {f_{z}^{t}\mspace{14mu}{\overset{\_}{f}}_{\overset{\_}{z}}^{t}} \right\rfloor = {{D\begin{bmatrix}\varphi^{t} & \psi^{t}\end{bmatrix}} = {D\begin{bmatrix}{\varphi_{k} + {t\;\Delta_{\varphi_{k}}}} & {\psi_{k} + {t\;\Delta_{\psi_{k}}}}\end{bmatrix}}}} \\{= {{D\begin{bmatrix}\varphi_{k} & \psi_{k}\end{bmatrix}} + {{tD}\begin{bmatrix}\Delta_{\varphi_{k}} & \Delta_{\psi_{k}}\end{bmatrix}}}}\end{matrix} & \left( {G{.11}} \right)\end{matrix}$where φ_(k), ψ_(k) are the current state of the (complex) variables andΔφ_(k) Δψk is the Newton search direction. The term on the left wasalready computed in (G.6). Hence, one only needs to compute the productD[Δφ_(k) Δψ_(k)] once, and for each step t, (G.11) boils down to ascalar-matrix multiplication followed by a (thin) matrix summation. Bothoperations are very efficient. The energy is then evaluated in adedicated kernel that accepts the matrix └ƒ_(z) ^(t)

┘ as input. Each thread handles a row that contains two complex numbers:ƒ_(z) ^(t),

, evaluate two real scalars: |ƒ_(z) ^(t)|²,

, and use it to compute the scalar E using Table 1. Finally, the kernelperforms a parallel reduction summation step to obtain E^(f) (Equation(31)).

In order to avoid CPU-GPU communication overhead, this step is carriedout on the GPU by using the dynamic parallelism feature which allowskernels to invoke other kernels, therefore allows the GPU to execute theline search autonomously.

The iterations are terminated when t∥Δφ_(k) Δψ_(k)∥ is smaller thanϵ=1e−12 (Algorithm 1, Line 9). The algorithm always converges.

In contrast, for some inputs, the SLIM and AQP methods struggle toconverge to machine precision error even after a thousand iterations.

Upon termination, the map itself is evaluated. This is done as explainedin connection with Equation (G.2). Since this product is also performedon the GPU, the result can be used by the GPU directly for rendering,avoiding expensive CPU-GPU transfer times.

For evaluating the performance of the inventive algorithm, the followingdefault parameters for all experiments were used, unless statedotherwise. The number of samples is |G|=10⁴, and |H|=0.1 |G|. The twoenergy terms E^(f) and E_(p2p) ^(ƒ) are balanced using λ=10⁵. Whencomparing to mesh-based methods, a mesh with 10, 000 vertices is used.User prescribed P2P target positions q_(i) are visualized as cyan disks.The images of the P2P under the map f(p_(i)) are visualized as smallerblack dots. A black dot which is centered inside a cyan disk indicatesP2P satisfaction. The computer runs Windows 10 with Intel Xeon (6 cores)CPU E5-2643 v4 3.40 Ghz, 32 GB, and has an NVIDIA TitanX graphics card.

FIG. 9 shows locally injective isometric deformation of Rex, adoubly-connected planar domain (left). Note the hole between the legs.The harmonic subspace has 392 DOF. The inventive solver converges after6 iterations to the result depicted, which is computed in 0.072 s whenimplemented on the GPU, and 0.375 s on the CPU. An unmodified Newtonsolver halts prematurely after 9 iterations at a high energyconfiguration due to a non-positive definite Hessian. The Newton-Eigensolver performs a full eigen-factorization Hessian modification at eachiteration. It converges after 32 iterations and takes 2.467 s which is×34 times slower than our GPU solver. The L-BFGS solver performs 1, 818iterations and is not competitive in terms of runtime.

FIG. 10 shows a comparison with BDHM using the same harmonic subspace.User distortion bounds for BDHM are set to: (k, σ₁, σ₂)=(0.8, 5, 0.3).The resulted map is bounded distortion with bounds (0.50, 4.3, 0.31).Without explicit distortion bounds, the inventive method minimizes thesymmetric Dirichlet energy producing a map with bounds (0.32, 2.25,0.43). The inventive GPU solver converges after 8 iterations with totalruntime of 0.0297 s. In comparison, BDHM requires 28 iterations ofsecond order cone programming, with total runtime 3.71 s, which is ×125slower.

FIG. 11 shows the effect of changing the number of Hessian samples |H|on the inventive algorithm when applied to the depicted deformation ofthe deer. Similar behavior was observed for other models anddeformations. With 200 samples or less, one sees an increase in thenumber of iterations needed for convergence (left plot). However, with400 samples, the convergence rate is identical to that of 50,000 (notethe three overlapping curves). On the right plot we see that with 50,000samples, the algorithm becomes 3 times slower due to the longer time ittakes to assemble the Hessian, yet, 400 samples are sufficient, in whichcase the algorithm converges after 30 ms.

FIG. 12 shows how the runtime of the inventive GPU-based modified-Newtonsolver is spread across the 4 main stages that compose a singleiteration. Each group of bins corresponds to a different model withvariable amount of degrees of freedom (variables). For all models,|G|=10,000, |H|=1,000. FIG. 12 demonstrates that the inventive solver iswell balanced with respect to the four different stages that compose asingle Newton iteration: gradient evaluation, Hessianmodification-and-assembly, Cholesky linear solve, and line search.

FIG. 13 (Table 2) shows a comparison of runtime (ms) of a singleiteration. SLIM and AQP use a mesh with 10, 000 vertices, while theother listed solvers operate within our harmonic subspace with thespecified DOF. Table 2 compares the runtime in milliseconds of a singleiteration of various algorithms minimizing the symmetric Dirichletenergy. Among all methods, AQP iterations are the fastest (followed byour GPU iterations), albeit we observe that AQP requires significantlylarger number of iterations compared to our GPU/CPU Newton solver.Moreover, the initialization required for AQP further slows down theruntime (right most column). For comparisons with AQP we run 5iterations of ARAP local/global followed by the projection method ofKovalsky et al, 2015 and 2016 (Kovalsky Shahar Z., Aigerman Noam, BasriRonen, and Lipman Yaron. 2015. Largescale bounded distortion mappings.ACM Transactions on Graphics 34, 6, Article 191 (2015); Kovalsky ShaharZ., Galun Meirav, and Lipman Yaron. 2016. Accelerated Quadratic Proxyfor Geometric Optimization. ACM Transactions on Graphics 35, 4, Article134 (2016)).

Throughout the experiments, the inventors consistently observed that theconvergence rate of the inventive algorithm is significantly better thanthat of the first-order methods (L-BFGS, SLIM, AQP) which is attributedto its second-order nature. Furthermore, the convergence rate is alwaysbetter even when compared against the Newton-Eigen method which employsa full eigen-factorization [Nocodal and Wright 2006, Equation 3.40] ofthe Hessian at each iteration, followed by a truncation of negativeeigenvalues to 1e-10. This is demonstrated in FIG. 9 (6 vs. 32 iters)and FIG. 14 (15 vs. 35 iters) suggesting that the Hessian modificationaccording to the invention is not only faster to compute but also moreeffective in reducing the energy. Both the eigen-factorization and thesolution to the linear system are computed using [cuSOLVER 2017].

The full potential of the inventive deformation algorithm is leveragedwhen it runs on graphics hardware, yet our CPU parallel Matlab version(which is approximately 5 times slower on average) is still much fasterthan the alternatives we compare to.

FIGS. 14 and 15 show a visual comparison of our algorithm againststate-of the-art methods minimizing E_(iso) ^(ƒ)+λE_(p2p) ^(ƒ). For eachmethod, small images of the maps obtained at intermediate stages of thealgorithm are included to show that for the competing methods, manyiterations are needed to get to a state where the progress is notvisually apparent. The inventive method is between ×25 to ×58 fasterthan these competing methods.

More particularly, the top three rows use the harmonic subspace, whileAQP and SLIM are general mesh-based methods. The total number ofiterations needed for each method as well as the total runtime isreported on the right of each row. The graphs depict the decrease inenergy as a function of iteration and as a function of time, separatingharmonic from mesh-based methods as they are scaled differently. The“jump” of 0.9 s in AQP's graph (first iteration) is due to theinitialization. Note that due to the large weight X, the energy at thefirst few iterations is very large and E_(p2p) ^(ƒ) obscures thebehavior of E_(iso) ^(ƒ). The zoom-in graphs show how the energydecreases more clearly. Our method converges in 15 iterations. Newton'smethod halts after 15 iterations due to non-positive definite Hessian.Newton-Eigen performs 35 iterations and is ×25 slower than our GPUsolver. AQP and SLIM require 400 and 120 iterations and are ×39 and ×58slower respectively.

Qualitatively, the results obtained with the mesh based methods tend toconcentrate the distortion near the positional constraints. Refining themesh does not seem to alleviate this issue completely while it increasesthe computation times even further. A mesh resolution of 10⁴ verticesprovides a reasonable balance between quality and speed and use it toevaluate AQP and SLIM. The inventive algorithm on the other hand israther indifferent to the mesh resolution and using a mesh with 10⁵vertices, leads to extremely smooth results with no degradation inspeed.

FIG. 16 shows additional results minimizing E_(iso) ^(ƒ) rendered athigh-resolution. Some of these deformations are intentionally radicaland are used to stress test the method and to demonstrate robustness.The inventive algorithm quickly converges on these extreme deformationsand certifies the maps as locally injective.

We claim:
 1. A computer-implemented method for interactive shapedeformation, the method comprising: obtaining an image; extracting ashape from the image, the shape being represented in a multiplyconnected planar domain comprising an exterior boundary curve (γ_(o))and one or more interior boundary curves (γ_(i)), the interior boundarycurves forming one or more holes of the shape; obtaining one or morepositional constraints; determining a planar shape deformation map,based on the one or more positional constraints, wherein the planarshape deformation map is determined by minimizing an isometricdistortion of the shape deformation, wherein the isometric distortioncorresponds to an integral over an isometric energy term; applying theshape deformation map to the shape to obtain a deformed shape; andoutputting the deformed shape, characterized in that the integralfurther corresponds to a boundary integral along the boundary curves(γ_(i)) of the shape and in that minimizing the isometric distortioncomprises checking whether the distortion is bounded on the boundarycurves (γ_(i)).
 2. The method of claim 1, wherein the shape deformationmap is locally injective.
 3. The method of claim 2, wherein the shapedeformation map is represented in an extended harmonic map space thatsupports shapes containing holes.
 4. The method of claim 2, wherein thestep of determining a deformation map is based on an analytic expressionfor the Hessian of isometric deformation energies.
 5. The method ofclaim 4, wherein the step of determining uses an approximation for theHessian of the isometric energy.
 6. The method of claim 5, wherein thestep of determining uses an analytic modification for the Hessian to bepositive definite.
 7. The method of claim 1, comprising the step ofestimating a flip-free condition for harmonic maps.
 8. The method ofclaim 1, wherein a Newton iteration is implemented on one or more GPUs.9. The method of claim 1, wherein the minimization is performed in areduced deformable subspace of harmonic maps.
 10. An article ofmanufacture comprising a non-transitory computer-readable medium havingprogram instructions stored thereon, the program instructions, operableon a device and, when executed by a processor on said device, cause saidprocessor to perform the method of claim
 1. 11. A device, includinghardware including at least one processor and at least one memory, thedevice programmed to perform the method of claim
 1. 12. The method ofclaim 1, wherein checking whether the distortion is bounded on theboundary curves (γ_(i)) comprises verifying that the shape deformationmap does not vanish on the boundary curves of the shape.
 13. The methodof claim 12, wherein verifying that the shape deformation map does notvanish comprises checking whether the map is orientation preserving overeach line segment of the boundary polygons.
 14. The method of claim 13,wherein checking whether the map is orientation preserving over eachline segment of the boundary polygons comprises checking whether theconditionσ₂(v _(i))+σ₂(v _(i+1))≥(L _(ƒ) _(z) +L _(ƒ) _(z) )l holds, wherein σ₂is the small singular value of the Jacobian, v_(i) and v_(i+1) areendpoints of the line segment, L are Lipschitz constants and l is alength of the segment.
 15. The method of claim 14, wherein verifyingthat the shape deformation map does not vanish further compriseschecking whether a turning angle of its holomorphic component changes onthe boundary curves.
 16. The method of claim 15, wherein checkingwhether a turning angle of its holomorphic component changes on theboundary curves comprises checking whether${{\sum\limits_{j = 0}^{N}{\sum\limits_{i}{{Arg}\frac{f_{z}\left( v_{i + 1}^{j} \right)}{f_{z}\left( v_{i}^{j} \right)}}}} = 0},$where the first summation is over polygonal loops representing holes ofthe shape, the second summation is over the segments in each loop, andwherein Arg is the principal branch of the complex argument function,f_(z) is the shape deformation map and v^(j) _(i) is the i-th point ofthe j-th loop.